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ABSTRACT 

An open question in the field of solar and stellar astrophysics is the source of heating that 
causes stellar coronae to reach temperatures of millions of degrees. One possibility is that the 
coronae are heated by a large number of small flares. On the Sun, flares with energies as low 
as those of microflares are distributed with energy as a power law of the form ^ cx E^°' with 
a ~ 1.8, and a appears to increase to values 2.2-2.7 for flares of lower energy. If the slope 
exceeds the critical value of 2, then in principle the entire coronal energy input may be ascribed 
to flares that are increasingly less energetic, but are more numerous. Previous analyses of flares 
in light-curves of active stars have shown that this index is generally > 2, though it may be as 
low as 1.6 when strong flares alone are considered. 

Here we investigate the contribution of very weak flares, covering the milliflare energy range, 
to the coronal luminosity of low-mass active stars. We analyze EUVB/DS events data from 
FKAqr, V1054Oph, and AD Leo and conclude that in all these cases the coronal emission is 
dominated by flares to such an extent that in some cases the entire emission may be ascribed to 
flare heating. We have developed a new method to directly model for the first time stochastically 
produced flare emission, including undetectable flares, and their effects on the observed photon 
arrival times. We find that apKAqr = 2.60 ±0.34, avi054Oph — 2.74 ±0.35, aADLco = 2.03 — 2.32, 
and that the flare component accounts for a large fraction (generally > 50%) of the total flux. 

Subject headings: X-rays: stars - stars: coronae - stars: flare - stars: late-type - methods: data analysis 
- methods: statistical 



1. Introduction 

The source of heating of solar and stellar coro- 
nae still eludes understanding even after decades 
of study (e.g., Schrijver et al. 1999). Despite sig- 
nificant evidence that magnetic activity is the 
prime driver for transferring energy into the 
corona, the mechanism by which this transfer oc- 
curs is not established in either the case of the Sun 
or other stars (see e.g., Rosner, Golub, & Vaiana 
1985, Narain & Ulmschneider 1996). Numerous 



heating mechanisms, such as acoustic wave dis- 
sipation (Stepien & Ulmschneider 1989), Alfven 
wave dissipation (Cheng et al. 1979, Narain & 
Ulmschneider 1990), magnetic reconnection phe- 
nomena (Parker 1988, Lu & Hamilton 1991) have 
been proposed, all of which might play some role 
in the overall heating. 

Recent work in the solar case has lent strong 
credence to the possibility of coronal heating being 
dominated by small-scale explosive events sugges- 
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tive of Parker's nanoRarc model, which is based on 
magnetic reconnections releasing energies ~ 10^'' 
erg event It has been well known that solar mi- 
croflarcs and milliflares"'^ are distributed in number 
as power laws of their energy output (Lin et al. 
1984, Hudson 1991), 



where E is the energy of the flare and fc is a con- 
stant. This relation has been verified and extended 
to lower energies by various authors, but despite 
the near universal acceptance of the form of the 
function in Equation 1 (e.g., Crosby, Aschwan- 
den, & Dennis 1993; see Kopp & Poletto 1993, 
Shimizu & Tsuneta 1997 for a different perspec- 
tive) , neither the index a nor the normalization k 
are well determined. For instance, a = 1.6 — 1.8 in 
the HXR to microflare energy range, and is vari- 
ously measured to lie in the range 1.8-2.9 at lower 
energies (Shimizu 1995 [a = 1.5 — 1.6], Porter, 
Fontenla, & Simnett 1995 [a = 2.3], Krucker & 
Benz 1998 [a = 2.3 - 2.6], ParncU & Jupp 2000 
[a = 2.0 - 2.6], Aschwanden et al. 2000 [a = 1.8], 
Winebarger et al. 2001 [a = 2.9 ± 0.1], Veronig 
et al. 2002 [a = 2.03 ± 0.09]). Recently Aschwan- 
den & ParncU (2002) have used scaling laws based 
on energy balance arguments to conclude that a 
must be ^ 1.6 on the Sun. The precise value of a 
is of considerable interest because if the power law 
is steep enough (a > 2), then in principle a mul- 
titude of small impulsive events would be suffi- 
cient to account for the energy output of the entire 
corona. 

Here we reconsider in particular an outstanding 
question in stellar X-ray astronomy, which is the 
nature of the apparently quiescent emission from 
active stars: does this emission actually arise from 
a superposition of a multitude of impulsive events 
(such as milliflares and microflares), or from truly 
quiescent plasma? Previous work based on de- 
tecting flares in EUV data (see e.g., Audard et al. 

^Because of the vast range of flare energies encountered, 
the energy ranges of the different flare types are not well 
defined. We adopt the convention (see e.g., Aschwanden 
et al. 2000) that milliflares cover the range E ~ io29-32 
ergs, microflajres E ~ j^g^®"^* ergs, and nanoflares E ~ 
2^q23-26 grgg \Yg consider all events down to the microflare 
regime to be 'normal' X-ray flare events, with similar origin, 
parameters, and effects, except for the differences in energy 
deposition. 



2000) suggests that flare contribution is indeed an 
important factor. Further, correlations of quies- 
cent X-ray flux with time-averaged U-band flare 
flux (Skumanich 1985, Doyle & Butler 1985) and 
the synchrotron radio luminosity (Giidel & Benz 
1993), together with the similar correlations found 
in the solar case (Benz & Giidel 1994) strongly 
suggest a link between the apparently quiescent 
emission and flares. In addition, spectroscopic 
evidence for high temperature plasma (T > 10^ 
K) during the quiescent phase (Butler et al. 1986, 
Kashyap et al. 1994, Drake 1996, Giidel et al. 1997, 
Giidel 1997) indicates that this quiescent emission 
could in fact be very similar to flare emission in 
origin. Thus, apparently quiescent coronae of ac- 
tive stars could be composed of a continuum of 
small unresolved flares, presumably distributed as 
power laws analogous to the Sun. This view is also 
supported by the double-peaked Differential Emis- 
sion Measures (DEMs) that result when an ensem- 
ble of flaring, hydro dynamically evolving loops are 
modeled on active solar analogs (Giidel et al. 1997, 
Giidel 1997). 

The possibility of stellar coronal heating due 
to small flares was considered by Ambruster, 
Sciortino, & Golub (1987) who searched for vari- 
ability in Einstein data of active stars and dis- 
cussed the contribution of low- level flaring to heat- 
ing stellar coronae. They concluded that while 
flaring must contribute at some level, the evidence 
does not justify extending the solar power-law dis- 
tributions to the stellar microflare case. Later 
studies of ensembles of strong stellar flares seen 
with EXOSAT and EUVE have shown these are 
distributed as power laws with index a = 1.6 — 1.8 
(Collura et al. 1988, Paflavicini et al. 1990, Os- 
ten & Brown 1999), thus ruling out low-intensity 
flares as a significant contributor to the heat- 
ing budget. In contrast, using a more sensitive 
method to detect fainter flares (sec Crawford et al. 
1970), Robinson et al. (1995, 1999, 2001) flnd that 
for stellar chromospheric and transition region 
events observed with the high-spcc;d photome- 
ter and the imaging spectrograph on the HST, 
a;~ 1.76 — 2. 17 in the chromosphere of the active 
dMc star CNLeo; a — 2.25 ± 0.1 in the chromo- 
sphere of the dMe star YZ CMi; and a ~ 2.2 - 2.8 
in the transition region of the dMOe flare star 
AUMic. (Note however that chromospheric and 
transition-region flare distributions have no known 
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direct correspondence with the coronal case.) Ap- 
plying a similar method to EUVE/DS data, and 
also correcting for overlaps in flares, Audard et 
al. (1999) find that for solar analogs EKDra and 
47Cas, a « 2.2 ± 0.2. This analysis was further 
extended by Audard et al. (2000) to a larger sam- 
ple of cool stars, and they find that a ranging from 
1.5 to 2.6, with the majority of the measurements 
having a > 2. Similar results are obtained for 
AD Leo (Giidel et al. 2001,2002). 

Note that the above studies are limited to rel- 
atively large flares {E > 10'^^ ergs) because of in- 
strument sensitivity, and also because the more 
numerous weaker flares are harder to detect in the 
presence of "contamination" by other weak flares. 
Thus, the low-energy end of the flare distribution 
is subject to large uncertainties. 

We have developed a new method to model 
the undetectable stellar flares and thus derive esti- 
mates of flare indices covering the milliflare regime 
as well as directly estimating the flare contribu- 
tion to the observed flux. We apply this method 
to active low-mass stars FK Aqr, V1054 Oph, and 
AD Leo. The datasets used are described in §2. 
The analysis method is detailed in §3 (a glossary 
of the terms used is given in Appendix A). The 
results are given in §4, and are summarized in §5. 



observed. Thus it is possible that flare heating 
could be a signiflcant component of coronal emis- 
sion for this star. Indeed, the EUVE/BS light 
curve shows evidence of a number of flare events 
(Figure 1) as well as large stochastic variability in 
the apparently quiescent emission. 

Table 1: FKAqr 



Other Names 


G1867A / HD 214479 


(RA, Dec)2ooo 


(22:38:45.56, -20:37:16.1) 


Spectral Type " 


dM2e/dM3e 


Period " 


4.08 days 


Distance ^ 


8.64 pc 


m.v 


9'".06 


B-V ^ 


1.47 


Lx " 


1.3 X 10^9 ergs s"^ 


EUVE/DS 


1997-oct-17 to 1997-oct-24 




(130.7 ks) 


Count rate 


0.36 ± 0.033 ct s-i 


Background 


~ 0.023 ct s-i 



"as in Strassmcicr ct al. (1993) 



L otrassmcicr ai. i^iyyo; 
'from the Hipparcos catalog (Perryman et al. 1997) 

1 the Einstein energy band (0.1-4.5 keV) (Dempsey et al. 
39.S1 



'•in 
1993) 



2. Data 

Here we analyze data from the Extreme Ultra- 
violet Explorer satellite's Deep Survey photome- 
ter (EUVE/BS) of 3 active low-mass stars.^ These 
stars are known to have significant flare activity, 
and do not undergo eclipses, and so are amenable 
to straightforward modeling: 

FK Aqr is a BY Dra type, spectroscopic, double- 
lined, non-eclipsing, low-mass, active binary (Ta- 
ble 1). Its flare energy output has remained steady 
over long intervals (8 years; Byrne et al. 1990), and 
optical modulation due to spot activity has been 



FK Aqr 



■^The intrinsic EUVE time resolution is 8 ms, and this is ad- 
equate to resolve the sources even at the maximum count 
rate seen in our observations (3 ct s~^). The Ei7VE/DS 
covers a useful spectral range of 52-246 A, with a peak ef- 
fective area of 28 cm~^ at 91 A. This wavelength range in- 
cludes many lines from highly ionized Fe XVIII to Fe XXIII 
normally found in the coronae of active stars, in addition 
to bound-free and free-free continua: plasma temperatures 
from Ri 1 to 30 MK are thus accessible for observation. 




200 300 400 500 
Time [ks] +927780.35113 
<500.0> 



Fig. 1.— EUVE/DS light curve of FKAqr. The 
light-curve corrected for instrumental effects is 
shown at a bin size of 500 s. The vertical bars 
denote the 1-sigma error on the count rate. Note 
that there are many obvious flares visible in the 
light curve, in addition to a base emission which 
also is highly variable. 
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V1054 Oph is a low-mass astromctric binary 
(Table 2). Previous X-ray observations with 
ROSAT have confirmed the existence of almost 
continuous flaring (Kcllct & Tsikoudi 1999) and 
the existence of high-temperature plasma that is 
responsible for most of the intensity variations 
(Giampapa ct al. 1996). EUVE data show a large 
number of relatively weak flare events that blend 
into the variations in the quiescent emission (Fig- 
ure 2). Such a datasct is very difficult to analyze 
by the traditional means of detecting and counting 
flares, but poses no difficulty to a direct modeling 
approach as is described below. 

Tabic 2: V1054 Oph 



VI 054 Oph 



Other Names 
(RA, Dec) 2000 
Spectral Type " 
Distance " 
mv " 
B-V " 

T b 

EUVE/DS 

Coimt rate 
Background 



Wolf 630 / G1644 / HD 152751 

(16:55:28.76, -08:20:10.8) 

M3Ve 

5.73 pc 

9^.02 

1.553 

5.6 X 10^^ ergs 

1994-jul-30 to 1994-aug-08 

(127.9 ks) 

0.09 ±0.014 ct s-i 

~ 0.015 ct s-i 



"from the Hipparcos catalog (Perryman et al. 1997) 
''in the Einstein energy band (0.1-4.5 keV) (Dempsey et al. 
1993) 



AD Leo is a well studied low- mass single flare 
star (Table 3) with a high flare rate. A long 
duration exposure was obtained by Giidel et al. 
(200L2002) that shows many large flares (Fig- 
ure 3). The data were obtained in 6 segments, 
and since the first segment could not be optimally 
reduced, and the last segment suffered from a high 
background, we have ignored them in this analysis 
and have concentrated on the 2"'* — 5*^* segments. 
In particular, we have carried out the analysis with 
the data grouped into two sets, segments 2 and 
3 forming one set and segments 4 and 5 forming 
the other. The light curves show slightly differ- 
ent characters in the two parts, with the former 
part dominated by large flares while the latter part 
shows smaller identifiable flares (Figure 3). 




)0 400 
Time [ks] +826330.28453 
<500.0> 



Fig. 2. — Same as Figure 1, but for V1054Oph. 
In this light curve it is difficult to cleanly distin- 
guish a flare event from the underlying continuous 
emission. 



Table 3: AD Leo 



Other Names 
(RA, Dec)2ooo 
Spectral Type " 
Distance " 
mv " 
B-V " 
r b 

El/VE/DS 



Count rate 
Background 



GJ 388 / SAO 81292 

(10:19:38.04, -M9:52:14.2) 

M3V 

4.9 pc 

9™. 43 

1.54 

8.91 X 10^8 ergs s"! 
1999-apr-05 to 1999-apr-14 
(258.6 ks) 

1999-apr-17 to 1999-may-04 
(347.7 ks) 

0.47 ± 0.03 ct s-i [apr05-aprl4] 
0.31 ± 0.03 ct s-i [aprl7-may04] 
~ 0.028-0.031 ct s-i 



"from Audard et al. (2000) 

''in the Einstein energy band (0.1-4.5 keV) (Dempsey et al. 
1993) 



In our analyses (see §3 below), we use the pho- 
ton arrival times directly, and apply the dead-time 
and Primbsching corrections^ for the particular 



^Priinbsching refers to the photons lost due to telemetry 
bandwidth, and is measured by the ratio of the total counts 
incident in a quadrant of the detector (as determined on 
the spacecraft; the summed counts from all the quadrants is 
used to determine the instrument deadtime) to the number 
of events telemetered to the ground (see the EUVE Data 
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AD Leo 



3.0 
2.5 
2.0 
" 1.5 

5, 1.0 



0.5 |f% 1ii>Wil)a,^«i||>i^^ 



400 60 
Time [ks] +973990.04492 
<500.0> 
AD Leo 



2.0 
„ '-5 
'" '-0 
" 0.5 

0.0 



-Xx X 



400 600 800 1000 
Time [ks] +975037.99421 
<500.0> 



Fig. 3. — Same as Figure 1, but for AD Leo; upper 

panel shows segments 2 & 3 (Apr 5 - Apr 14), and 
the fower panel shows segments 4 & 5 (Apr 17 - 
May 4). Many large flares arc visible (especially 
in segment 2, which is dominated by a large flare), 
as well as significant variability in the apparently 
quiescent emission. 

observation to the models (sec §3.1) over good 
time intervals (GTIs) defined to exclude SAA pas- 
sages and Earth blockages. That is, we process 
model light curves to produce a simulated set of 
photons and apply time windowing and statisti- 
cal censoring (i.e., the instrument response in the 
time domain) to generate an event list that may be 
directly compared to the observed events. For the 
sake of brevity, these corrections will henceforth 
be referred to as "instrumental" corrections. The 
source photons are collected over a circle of radius 
4" surrounding the source. In all cases, < 10% of 
the events are estimated to be due to backgroTind 
photons, except in the case of V1054Oph where it 
is estimated to be 17%. 

3. Analysis 

Previous attempts to determine the values of 
the parameters in Equation 1 have concentrated 
on first detecting flares in the binned light-curve 
and then fitting power law expressions to the de- 
tected numbers. In contrast, we assume the reality 
of a power law distribution and set up a model to 
compare with the data. This model is described 



Products Guide for a detailed description of its origin and 
correction). 



in §3.1, and the manner in which the model pa- 
rameters are derived is described in §3.2. The ap- 
plicability of the method, including verification, 
assumptions, advantages and disadvantages, are 
discussed in §3.3. 

3.1. The Model 

Because flares generally occur randomly (see 
§3.3.2), we cannot directly model the light curves, 
as it is not possible to "deconvolve" a complex 
light curve by specifying the location of each flare 
in a model (see, e.g.. Figure 2). Instead, rather 
than match the flare locations in detail, we carry 
OTit a fitting process wherein only the number and 
intensity distributions of a set of model fiares are 
matched with the data. This is accomplished by 
comparing the distributions of photon arrival-time 
differences. The assumptions made in defining the 
model described below are discussed in detail in 
§3.3.2. 

We first generate a set of photon arrival times 
by simulation from a 3-parameter model 



M = {a,rF,rc} 



(2a) 



where a is the index of the power law as in Equa- 
tion 1, rp is the average count rate due to flares, 
and rc is a constant component which is also ex- 
pected to fully account for the background (see 
§3.3.2). It is also possible to use the average to- 
tal flux vt = rp + rc {rp < tt) as the defining 
parameter instead of rc, 



M' = {a,rp,rT}, 



(2b) 



with equivalent results. The counts at time t in an 
interval dt, C(t) ~ Poisscm[T(t)dt] are Poisson- 
distributed according to the instantaneous rate 
r{t). The rate r{t) may be described as due to 
the sum of model counts due to a flare compo- 
nent f{t) and a non-flare component rc{t), and 
a correction factor 4){t) that takes into account 
Primbsch, dead-time, and GTIs,'', i.e., 

C{t) ~ (l){t) Poiss(m[rc{t) dt + f{t) dt ] , (3a) 

where rc{t) is taken to be constant unless stated 
otherwise. The flare component f{t) may in turn 



^ Note that < </>(i) < 1 determines the probability that 
any photons are collected at time t, and in particular is 
identically outside the GTIs. 
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be written as a superposition of numerous individ- 
ual flares, i.e., 



Nf 



f{t) = Y,^{t-t,)F,e 



it-t,)/r 



(3b) 



where r is the flare decay timescale,^ Nf are the 
total number of flares, Fj are the peak energy in- 
tensities of individual flares (the counts due to 
each flare, c = tFj are sampled from the distri- 
bution represented by Equation 1), and 



e(x) 



X < 
a; > 



is a step function to represent flare onset. 

Note that not only will the placement tj and 
peak intensity Fj of the flares vary for each sim- 
ulation, but so will the total number of flares Nf. 
Within the bounds of Poisson statistics, we expect 
that for any given simulation, the total energy due 
to all the simulated flares is determined by the av- 
erage expected count rate due to flares and the 
total duration of the observation, AT, i.e.. 



/ dtf{t) = J2FjT = rFAT. 



(3c) 



Note that Equation 1 is written as a function of 
the energy deposited by a flare E, but assuming 
that the observed counts due to this flare c is pro- 
portional to E (see discussion in §3.3.2), i.e., 



dN dN 
dE°^~d^ 



'dc. 



(4) 



we can use the model parameter rp to fix the nor- 
malization K of the power law. By equating the 
total counts due to the flare component with the 
counts expected from the power law distribution, 
we get 



/ dec"'' 



dc 



AT 



(5) 



The upper limit in the integration is deflned by 
requiring that all the observed counts be due to a 
single model flare {cmax = max{Fj t} = rp AT), 

3 



^We assume r to be fixed for purposes of simplicity. See 
§3.3.2 for a discussion of cases when t may vary. 



that is, no flare model may produce a flare with 
more counts than are observed. The lower limit in 
the integration is deflned by requiring that each 
flare be assigned at least 2 counts (i.e., Cmin — 2; 
this is so that an arrival time difference may be 
determined even in the extreme case where the 
model may have just one weak flare see §3.2 be- 
low). Thus, carrying out the integral in Equation 5 
and rearranging the terms, 



«l(a#2) 

K|(a=2) 



{2-a)rpAT 



{rpATy-<^ - 

rpAT 
ln{rpAT/2) 



22- 



(6a) 
(6b) 



For a > 2, this implies that if « rr, then the 
adopted lower limit is very close to the theoretical 
lower limit to the extent of the power law distribu- 
tion in order for it to account for all the observed 
counts (see Table 4). 

In order to obtain the best values of the pa- 
rameters that describe the data, and a confidence 
range on these parameters, we carry out a forward- 
fitting procedure based on a Bayesian formalism 
(see e.g., Loredo 1990 for a tutorial on the founda- 
tions of Bayesian probability theory): we compute 
the probability distribution of the model parame- 
ters given the data, which is a composite of what- 
ever prior information we may have on the model 
parameter values and the likelihood of realizing 
the observed data for specified parameter values. 
That is, we derive the joint posterior probability 
p{M.\D, I) of the model parameters conditional on 
the data,^ 



p{M\D,I) oc 



p{a\I)p{rp\I)pirc\I) 
xp{D\M,I), (7a) 



where D represents the data, and / represents as- 
sumptions necessary to solve the problem, includ- 
ing the effects of instrument characteristics. The 
first 3 factors on the right hand side of the equa- 
tion are the a priori probability distribution func- 



^The expression p(x) represents the probability that the log- 
ical statement "x" is true. In particular, conditional state- 
ments are written as i.e., p{A\B) represents the 
probability that statement "A" is true given that state- 
ment "_B" is true. These statements may be generalized 
to include models and parameter values; as an illustrative 
example, p{{a = 2.1,rp = 0.1, = 0.3}|D,/) = 0.1 reads 
as "the probability is 0.1 that a = 2.1, rp = 0.1, and rx = 
0.3 given the data D and supporting information I." 
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tions of the model parameters, and the last fac- 
tor, p{D\'M.,I) is the likelihood of obtaining the 
observed data given the model parameters. Com- 
puting the posterior probability distribution con- 
stitutes a complete solution to the inference prob- 
lem for the specified model M. Note that the above 
is a simplified form of Baycs' Theorem wherein the 
normalization factor p{D\I) usually present on the 
right hand side of the equation is ignored (e.g., 
Kashyap & Drake 1998, van Dyk et al. 2001). If 
Tt is used instead of rc, the joint probability dis- 
tribution takes the form 

p{M'\D,I) cx p{a\I)p{rF\rT,I)p{rT\I) 

xp{D\M',I), (7b) 

with p{rF\rT,I) = for rp > Tt- In the fol- 
lowing, we make no distinction between M and 
M'. The Bayesian formalism allows us to deter- 
mine the probability distributions of each of the 
parameters by marginalizing, i.e., integrating over 
the other parameters. Thus we can write for the 
probability distribution of a alone, 

p{a\D,I)= I drp f drc p{M\D,I) (8) 

and similarly for the other parameters (cf. Equa- 
tion 11). 

3.2. The Algorithm 

The problem facing the modeling process is il- 
lustrated in Figure 4, where the light curve from 
the observation of FKAqr is compared with se- 
lected simulated model light curves.'' It is easy to 
recognize that the model light-curve for a = 2.5 
is the most similar in character (i.e., in the num- 
ber and strength of discernible flares) to the data 
light-curve, but clearly the locations of the flares 
do not match. Normal fitting methods that rely on 
matching the expected model counts in a bin with 
the observed counts would fail on a problem such 
as this. We thus seek to employ a method which 
compares the interesting information between the 
data and the model without being misled by the 
obvious, though uninteresting, differences. 

One method that would satisfy our require- 
ments of simultaneously ignoring flare locations 



FK Aqr: EUVE/DS 



'^Wc emphasize that these Ught curves arc shown only for 
purposes of illustration, and that the analysis does not re- 
quire binning the observed events. 
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Fig. 4. — Illustrative comparison of data and 
model light curves. The light-curve from an obser- 
vation of FK Aqr is plotted in the topmost panel 
at a bin size of 500 s. No instrumental corrections 
have been applied to the data. In the lower pan- 
els, light-curves derived from simulated events for 
various values of the index of the power law dis- 
tribution, a = 1.8,2.1,2.5,3.0 are shown at the 
same binning, and after applying the appropriate 
Primbsch and dead-time corrections and GTI fil- 
tering. The model curves are computed assuming 
rc = in order to simplify the comparisons. Note 
the larger dynamic range in the light curves for 
lower a. It is apparent that the light curve for 
a = 2.5 is most "similar" to the observed light 
curve. 

and yet be sensitive to flare numbers and inten- 
sities is to compare the distributions of photon 
arrival-time differences g{6t) (where 6t is the inter- 
val between consecutive counts) between the data 
and the model.® In the absence of any variability, 
i.e., if the light curve is flat with expected rate r, 
the resulting set of St are distributed as an expo- 



Other methods such as computing the fractal length of the 
events, multi-scale analyses of light curves, etc., may also be 
applied (see e.g. , Vlahos et al. 1995) . A detailed comparison 
of the benefits of one method versus another is beyond the 
scope of this article, but note that the results we derive 
here axe robust within the regime of applicability of the 
adopted method (see §3.3). 
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Fig. 5. — Comparison of distributions of arrival 
time differences between the data and model. 
Note that the x-axis is in log-scale. The upper 
plot shows the full distributions g{St) for the data 
(stepped curve), a = 1.8 (dotted curve), a = 2.1 
(dashed curve), a = 2.5 (solid curve), and a = 3.0 
(dash-dotted curve). The models used are the 
same as in Figure 4. It is apparent that the first 
two models are bad fits to the data. A more de- 
tailed analysis is required to select between the 
last two models. In order to highlight the differ- 
ences in the distributions near their peaks, these 
differences are shown in the lower plot, where the 
differences between the model and the data distri- 
butions are plotted for a = 2.1 (dashed), a = 2.5 
(solid), and a = 3.0 (dash-dotted). A comparison 
of the values indicates that a = 2.5 is a signifi- 
cantly better fit to the data than a = 3 (x^ ~ 465 
and 490 respectively). 



nential, 



g{5t) are 



-r St 



(9a) 



In the presence of variability in the expected rate, 
the observed distribution would be a superposition 
of many such distributions: if the fraction of time 
that a source spends at rate r, is given as pi, then 



g{5t) (x^piriC' 



(9b) 



ent magnitudes of variability. In particular, events 
generated from a model with low a (e.g., a = 1.8, 
where the light curve would be dominated by a 
few very large flares) would result in a distribution 
g{St) that is skewed to smaller values of St, while 
those from a model with large a (e.g., a = 3.0, 
where the light curve would be composed of a 
large number of very small flares that overlap each 
other so finely that it would not be possible to dis- 
tinguish it from a source with constant intensity) 
would approach the limiting case of Equation 9a 
above. 

In Figure 5 we show the comparison between 
the distributions of arrival-time differences derived 
from the same datasets in Figure 4. The differ- 
ent curves may be compared using any of a num- 
ber of statistical methods such as computing the 
X^, applying the Kolmogorov-Smirnoff test, etc. 
This approach succeeds in providing us with an 
objective measure of the "similarity" between the 
datasets; indeed of the 4 models considered in Fig- 
ures 4 and 5, the one with a = 2.5 has the smallest 
X^ when compared with the data. Note that Fig- 
ure 5 also illustrates a fundamental limitation of 
this method, viz., the method loses sensitivity for 
larger values of a, which may be indistinguishable 
from sources with a constant intensity (see §3.3.1). 

Because the model (Equation 3a) is stochastic, 
we use Monte-Carlo simulations to generate many 
realizations for each set of model parameters. The 
simulated distributions of arrival-time differences 
gMODEhiSt) are compared with the corresponding 
distribution derived from the data go at a (St) over 
a parameter grid. The likelihood is computed as 
the probability density of obtaining the observed 
X^ value for N degrees of freedom (see Eadie et al. 
1971, their Equation 4.22): 



p{D\M) 



1 ('x^ 

2\2J 



r(f) 



(10) 



It is thus possible to distinguish between differ- 



The a priori probability distributions for the pa- 
rameters (Equation 7a) are taken to be non- 
informative, and thus flat, over the range of pa- 
rameter values defining the grid. The basic steps 
in the algorithm we follow are outlined below: 

1. From the data, derive the distribution of 
photon arrival-time differences, excluding 
the gaps in the data due to breaks in the 
GTIs. 
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2. For the specified values of model parameters 
M (Equation 2), obtain a realization of the 
photon event-list over the duration of the ob- 
servation. This is done by first simulating a 
light curve incorporating all the flare events 
(Equation 3b), added to the base emission, 
and then deriving photon arrival times based 
on the instantaneous count rates. 

3. Because we are to compare EUVE/DS event 
list data with simulated data, instrument ef- 
fects must be taken into account. This is 
encoded in the factor (j){t) (see Footnote 4). 
We apply Primbsching and dead-time cor- 
rections by discarding photons with proba- 
bility 1 — 4>{t) (i.e., sample a random number 
z from a uniform distribution over the inter- 
val [0, 1], and discard the photon at time t 
if 2; > 4'{t)). The retained set of events is 
identical in its instrumental characteristics 
to the observed data. 

4. Prom this set of simulated photon arrival 
times, derive the model distribution of 
arrival-time differences over a binning that 
maximizes the reduced x^-^ 

5. Compute the likelihood as in Equation 10 

and the a posteriori probability at the spec- 
ified grid point as in Equation 7a. Note that 
the number thus obtained is not normalized, 
and so must not be used to compare, for ex- 
ample, the relative probabilities of different 
types of models. 

In practice, the above algorithm must be en- 
hanced by some additional steps. For instance, 
the likelihood may be artificially reduced if a large 
model flare is fortuitously located coincident with 
large dead-time. We therefore shift the Primb- 
sching and dead-time corrections by a random in- 
terval and recompute the arrival-times. This is 
mathematically equivalent to shifting the simu- 
lated events, but is done in this fashion because 



We do not know the optimal binning for g{St) a priori. The 
binning must be chosen such that differences between the 
two distributions being compared are highUghted to best 
advantage, at a scale that is determined by the datasets 
themselves. Because the total number of photons in the 
datasets being compared are approximately the same, very 
coarse binnings and very fine binnings both produce low 
values of the reduced x^i ''•nd we adopt as the optimal 
binning that which maximizes this value. 



of the lower computational cost. Typically 3 
such shifts are carried out for each simulated light 
curve, and the best comparison is chosen. In ad- 
dition, in order to derive a robust estimator, we 
perform f» 15 — 20 simulations for each M (or M') 
and adopt the median value of the set as the final 
value of ^(MII?, /). The grid of model parameters 
are chosen such that a is explored over the useful 
range of the algorithm (1.5. . .3.0; see §3.3.1); Tt in 
a range within Scr of the average count rate; and 
rc and rp ranges from « to the average count 
rate. There are typically ~ 15 grid points for each 
parameter. 

The probability distribution along any of the 
axes is then obtained by summing over the other 
axes and normalizing, e.g., 



p{a\D,I) 



EEp(m|A/) 

rp rc 



(11) 



The derived probability distributions may be sum- 
marized by their means and variances, e.g., 

E ap{a\D,I) 

a = 



var{a) = o? — {of' . 



(12a) 
(12b) 



3.3. Applicability 

As was demonstrated above (§3.2, Figures 4 & 5), 
the distributions of arrival-time differences g{dt) 
provide an objective means to compare event lists 
that are dominated by stochastic events. Here, 
we verify that the algorithm gives reasonable re- 
sults by creating simulated datasets and obtaining 
best-fit values for them (§3.3.1), then discuss the 
effects of some of the assumptions we have made 
in formulating the problem (§3.3.2), and detail 
the advantages and disadvantages of the adopted 
method (§3.3.3). 



Other well-known methods of summarizing probabil- 
ity distributions include noting the MAP (maximum a 
posteriori value; the mode of the distribution), as e.g., 
p(ay,.p\D, I) = max{p(a\D, I)}, or a range correspond- 

a 

ing to an integrated area under the curve, equivalent to 
some defined probability vr = J p(a\D,I), such that 



G [a. 



etc. Unless otherwise specified, we 



always report the mean values and la errors. 
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3.3.1. Verification 

Wc show here that the algorithm works as ex- 
pected by simulating well defined datasets and 
then "fitting" to them. We generated flare models 
with specific values of a (labeled ol^j^^j^) ranging 
from 1.5 to 3, and for the sake of simplicity, held 
Tp = 0.5 ct and rc = 0.1 ct s~^. Simulations 
were carried out over a time interval AT = 100 
ks and assuming a fixed flare decay timescale of 
T = 3 ks. These values were chosen as being typ- 
ical of EUVE/DS observations. The datasets to 
be fit to were chosen by first simulating 9 sepa- 
rate event lists and then choosing the one with the 
median number of counts, in order to avoid con- 
taminating the verification process with extremal 
datasets. The fits were carried out as described 
in §3.2, and the resulting best-fit values 
shown in Figure 6. The fitted values follow the 
true values closely, for all 3 parameters. 




1.5 2.0 2.5 3.0 

'"true 



Fig. 6. — Verifying the method. The best-fit val- 
ues a^jr^ obtained for simulated datasets gener- 
ated using a specific values a^^^^^, are shown as 
diamonds, and the 90% credible regions are shown 
as vertical bars. The line of equality (solid line) is 
also shown. 

The method works well over the range of inter- 
esting a, that is, spanning the value a — 2, and 
is capable of distinguishing between sources with 
a above or below this critical value (see §1). The 
sensitivity of the method naturally decreases as a 
increases to > 3 when the datasets become indis- 
tinguishable from that of a steady source. How- 
ever, this is not a hard limit and can be extended 
for datasets with larger average flare component 



intensity rp and longer observations AT. The 
reliability also decreases as a decreases to < 1.5 
since at these values of a the simulations are dom- 
inated by a very few but very large flares and are 
therefore subject to large fluctuations, and hence 
are not robust. A larger number of simulations at 
each grid point becomes necessary when a dataset 
exhibits smaller values of a. 

3.3.2. Assumptions 

We have made a number of simplifying assump- 
tions in our analysis, and these are discussed be- 
low, with particular attention to the effect they 
have towards the robustness of the results. In gen- 
eral, all our assumptions are conservative, in the 
sense that they all act to generate a best-fit a and 
rp that is smaller than the true a and rp. That 
is, the main results expounded here, that the flare 
distributions on active stars have a > 2, and that 
the apparently quiescent emission is dominated by 
flares, are not affected. 

Power laws: Energy release due to flare events 
in both solar and stellar environments has been 
well established to be highly intermittent and that 
the events are distributed as power laws spanning 
many decades in energy (e.g., Crosby et al. 1993, 
Giidel et al. 1997, Giidel 1997, Krucker fc Benz 
1998, Osten & Brown 1999, Audard et al. 2000, 
Aschwanden et al. 2000, Veronig et al. 2002). This 
suggests the absence of a characteristic scale for 
the intensity of a flare event, and is understood 
to arise from avalanche or SOC (self-organized 
critical) models (Lu & Hamilton 1991, Vlahos et 
al. 1995, Georgoiflis & Vlahos 1998, Krasnosel- 
skikh et al. 2001). Nevertheless, there is evidence 
from studies of solar flares that the index of the 
power law distribution does not remain the same in 
the microflare and nanoflare range (cf. Aschwan- 
den et al. 2000 [their Figure 10], Winebarger et al. 
2001). Studies of stellar flares (e.g., Ambruster et 
al. 1987) also suggest that more than one type of 
plasma instability may be present, and that the 
distribution may steepen or change at lower flare 
energies; for instance, compare a « 1.8 found by 
Osten & Brown (1999) with the generally larger 
values found by Audard et al. (2000) who include 
much weaker flares in their analysis (also see the 
results from AD Leo presented here in §4). Stud- 
ies of flare distributions arising from the transition 
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region also show a similar dichotomy (Robinson et 
al. 2001). Such changes in power law indices could 
arise from a variable driving mechanism (Geor- 
goulis & Vlahos 1998). It is therefore possible 
that the true distribution departs from exact self- 
similarity in some complicated manner. However, 
present data, especially in the case of stellar flares, 
are insufficient to detect these departures (Audard 
et al. 1999,2000; also Giidel et al. 2001,2002). Here 
wc assume that a single power law index is valid 
over at least 4 orders of magnitude (see Table 4) 
in flare energy. If the distribution does steepen 
for lower flare intensities, then note that first, the 
steepest parts, which approach the limiting case 
of constant emission (see Equation 9a), will con- 
tribute to enhancing the constant intensity rc (as 
described in §3.2). Hence the fitted rp will be a 
lower limit to the true value. Second, the fitted 
a will be a weighted average biased towards the 
high count rate, shallower distribution, and inas- 
much as a "true" value of a may be said to exist, 
it would be greater than the fitted value. 

Decay Timescales: We analyze the data by fit- 
ting the 3 parameters a,rp, and rc of the model 
(Equation 2). We assume that the flare decay 
timescale r (Equation 3b) is fixed (usually at 3 
ks, as suggested by the detectable fiares, and the 
radiative cooling timescales suggested by ROSAT 
observations; see e.g., Giampapa et al. 1996) and 
is the same for all flares. However, it is well 
known that t varies for individual flares on ac- 
tive stars (e.g., Pallavieini et al. 1990) and espe- 
cially so for RS CVn stars (Osten & Brown 1999). 
Note though that we confine ourselves to a spe- 
cific class of active stars low-mass main-sequence 
stars that flare frequently - and exclude RS CVn 
stars from our sample. There is evidence that 
flare duration scales with flare energy in various 
passbands (e.g., Crosby, Aschwanden, & Dennis 
1993; also, Vlahos et al. 1995 for avalanche mod- 
els, Temmer et al. 2001 for H^ flares, Georgoulis, 
Vilmer, & Crosby 2001 for deca-keV flares, Robin- 
son et al. 2001 for chromospheric and transition re- 
gion events, Veronig et al. for GOES SXR flares) 
such that more intense flares appear to last longer. 
Based on avalanche model simulations, Lu et al. 
(1993) flnd that on average, an event with de- 
cay timescale r corresponds to an energy release 
E oc T^-^^. (See also the discussion below about 



sympathetic flaring, which could affect measure- 
ments of flare durations.) But note that the decay 
timescales for soft X-ray flares (such as the ones 
wc arc concerned with) are primarily dependent 
on flare temperature and plasma density, and sec- 
ondary effects such as changes in the heating rate 
come into play only for very large flares where the 
heating timescales approach or surpass the radia- 
tive cooling timescales. In a model that incorpo- 
rates such variations of r, the fainter flares will 
have larger peak rates than in the regular model 
(in order to have the same total energy output), 
and the distribution of arrival-times g{St) will be 
skewed towards smaller dt. Thus, fltting to data 
that may be generated in this manner using a r 
flxed by the higher intensity flares would result in 
a^jr^ < a^nuE (s6e Figure 5) and therefore the flt- 
ted values would be lower limits to the true indices. 
Conversely, fltting a model where r decreases with 
flare energy to a dataset that is not generated in 
this manner will result in a^^^ > a^^^^^. (See 
also discussion by Giidel et al. 2001,2002.) In or- 
der to test the sensitivity of our datasets to these 
variations in modeling, we have carried out fits to 
the data using models with an energy-dependent 
timescale, r oc , in particular, P = j (Giidel 
et al. 2002, based on fits to EXOSAT flare decay 
timescales of Pallavieini et al. 1990). We flnd that 
the best-flt value of a increases by ~ 0.3 — 0.4 for 
the complex model, e.g., for segment 3 of AD Leo 
a changes from 2.2 to 2.7 . We thus conclude 
that using the simpler model (r = constant) is 
preferable in that we do not overestimate a, and 
further note that stronger dependences, e.g., that 
suggested by the theoretical study of Lu et al. 
(1993; T oc E^-^^) would result in even larger 
fltted values of a. In addition, we have also ex- 
plored the sensitivity of flts to the adopted value 
of T, by generating simulated data with small de- 
cay timescales t,,jj^j and then fltting to it a model 
with a larger decay timescale Tpj^^,. We find that 
as expected the best-fit a is smaller than the true 

value. Ori^lir >r 1 < Q^qta^It -^^ 

^This effect is large only for small values of a. For ex- 
ample, {a^jj,^ = 1.7;t^;^ = 1000} is fit by {a^^^ = 
1.5 ± 0.02; r^^j,^ = 3000}, whereas {2.1;500}sj„ is fit by 
{2.16 ± 0.17; .3000}^jy (i-e., is indistinguishable from the 
true value), etc. Note that, as expected, this exercise also 
shows that when the data are dominated by multiple over- 
lapping flare events, the results are insensitive to the pre- 
cise value of the decay timescale, and effects such as stellar 
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The reverse holds true for the opposite case, i.e., 
ctpr-rlfT <T 1 > ct«!r»^lT ! ih which case 
the best-fit values imply flare distributions that 
are steeper than they should be. Thus, the rela- 
tively large value we adopt for r (3 ks) provides 
a conservative estimate, and the fitted values may 
be considered lower limits to the true values. 

Rise Times: We model the flares as having an 
instantaneous rise and a slow exponential decay 
(Equation 3b). In general this is a good approxi- 
mation, since the rise times are short compared to 
the flare duration (e.g., Reale, Peres, & Orlando 
2001, Temmer et al. 2001). Further, our algorithm 
is insensitive to the exact sequence of the emis- 
sion intensities (i.e., the process of forming g{5t) 
destroys the sequential information in the event 
lists), and thus any discernible flare rise times will 
be indistinguishable in their effects from flares of 
small decay times. Thus, this approximation re- 
duces to the problem discussed above, that of a se- 
quence of flares that are "contaminated" by flares 
with smaller decay timescales, with similar effects 
on a and rp- 

Flare Waiting-times: Much work has been 
done to characterize the time between flares within 
individual active regions. If we consider each flare 
to be an independent event then the waiting time 
between flares is a Poisson process (Rosner & Va- 
iana 1978). However, actual observations of solar 
Hard X-ray flares show that the waiting-time dis- 
tribution (WTD) is a power law in intervals, with 
index ranging from —2.16 to —2.4 (Boffetta et 
al. 1999, Wheatland 2000, Lepreti, Carbone, & 
Veltri 2001), though Moon ct al. (2001) flnd that 
for strong solar flares (strength greater than CI) 
the waiting times are indeed well characterized by 
a Poisson distribution. The power law distribu- 
tion has been interpreted as due to sympathetic 
flaring (i.e., a cascade of small flares that follow 
a large flare, thus invalidating the assumption of 
event independence) by Wheatland, Sturrock, & 
McTiernan (1998), and adapted within SOC mod- 
els as a non-stationary random process (Norman 
et al. 2001). However, since we can only observe 
disk-integrated flux in stellar data (it is not pos- 
sible to monitor individual active regions), event 

rotation may be ignored. 



independence is a better approximation. We thus 
assume in our modeling that the stellar flare WTD 
is Poisson. In such a case, the model undercounts 
the number of flares separated by short intervals 
(Wheatland et al. 1998). These flares would gen- 
erally be associated with the stronger flares that 
set off a cascade, thus effectively increasing the 
decay timescale for large flares. As argued above, 
this causes ctpj^ < ot^^^j^. 

Energy Deposition: We have implicitly as- 
sumed that the observed counts track the energy 
deposited by the flares linearly (see Equations 1 
and 4). That is, we assume that the energy depo- 
sition process that causes the flare event (gener- 
ally considered to be magnetic reconnection - see 
Parker 1988) is, first, sufficiently energetic that a 
large fraction of the deposited energy goes towards 
thermal loading of the plasma and not into bulk 
motions (see Winebarger et al. 2001); second, that 
the resulting plasma temperatures for all flare en- 
ergies lie near 10^ K; and third, that the tempera- 
ture evolution of the plasma after the flare event is 
not drastic. These assumptions are supported by 
the emission measure analysis of several Yohkoh 
flares by Reale et al. (2001), who find that solar 
flares are dominated by emission at ^ 10 MK. 
Hydrodynamic modeling of an ensemble of flaring 
loops by Giidel et al. (1997) and Giidel (1997) 
also shows that the bulk of the flare DEM lies 
above 6 MK; indeed Giidel (1997) finds that the 
DEM is bimodal around 10 MK, attributable to 
the slightly smaller temperatures generated by the 
weaker flares (but which are still significantly hot- 
ter than the temperatures achieved by the quiet 
Sun). Thus, while it is reasonable to expect that 
flare emission would evolve from temperatures 
of ~ 20 MK to - 5 MK, and that flare events 
that are less energetic would heat the plasma to 
a lower temperature, in practice the effects of 
such variations are very little. Furthermore, be- 
cause we model the distribution in counts space, 
i.e., ^ rather than ^ directly, we invariably 
obtain distributions that are shallower than the 
true distributions^^ (see also extensive discussion 

'Suppose wc write the observed counts generated in a detec- 
tor Cabs due to emitted energy Etrue as a power law tiiat 
deviates from linearity by a small amount, c^bs ^trus^ 
(5 > 0. That is, as Etme decreases, the counts produced 
in the detector decrease at a faster than linear rate, as 



12 



in Giidcl ct al. 2002). In the EUVE/DS, changes 
in temperature of this magnitude causes the ob- 
served counts to vary by a factor of w 2. Fur- 
thermore, the range of flare energies we can model 
(see Table 4) are much larger than the relatively 
low-energy "explosive events" that are character- 
ized by large non-thermal velocities (Winebargcr 
et al. 2001 and references therein), and which may 
not contribute significantly to the plasma heating. 
That is, we assume that lower energy depositions 
that may heat the plasma to lower temperatures 
would not be detectable by the EUVE/DS in any 
case. Note that the assumed abundances will also 
affect this to some extent, but its effect is minimal 
for a broad-band instrument such as is used here. 
We therefore assume that an observed EUVE/DS 
count corresponds to a photon of average energy 
1.7 X 10"" ergs cm-^ ct"^ over the 59-250 A 
range produced by a plasma at 10^ K (PIMMS 
v2.5). 

Cut-offs: We have adopted upper and lower lim- 
its to the extent of the power law distribution (see 
Equation 5) that are based strictly on numerical 
expediency: the upper limit is set by the require- 
ment that the largest possible flare can produce 
no more than the observed number of counts, and 
the lower limit is defined by the minimum num- 
ber of counts needed to define an interval. How 
do these limits compare with theoretical expec- 
tations? First note that the solar flare distribu- 
tion has been explored to much lower energies, 
E - lO^'^ ergs (Winebarger ct al. 2001), than has 
been achieved for stellar observations. Recently 
Katsukawa & Tsuneta (2001) have estimated that 

would be expected to happen if the lower-energy flares re- 
sult in lower plasma temperatures, and the detector has 
a smaller response to these temperatures. This is a rea- 
sonable approximation for broad-band instruments such 
as the EUVE/DS, though it may not be applicable for 
small-passband detectors such as TRACE (see Aschwan- 
den & Charbonneau 2002). If j^f^ oc E~"l'-"' , then 

-£~~ °c , where cio^s = "'i^V"'^ » ^'^'i always 

smaller than atrue for <5 > and atrue > 1- Conversely, 
if 5 < 0, as may happen for extremely high flare energies 
and temperatures that lie above the best response of the 
detector, the observed distribution will be steeper than the 
true distribution, a^ts > citrue', however this case has little 
effect on our analysis results because of the wide temper- 
ature response of the EUVE/DS, and any flares that fall 
outside its range would be too few in number to affect the 
results. 



the probable energy of Parker-type nanoflares is 
< 10^^ ergs. In contrast, we find that if the ob- 
served EWE emission is assumed to originate en- 
tirely in flares, then it is sufficient to extend the 
power law distributions to energies E ~ X0^^~^^ 
ergs (see Table 4). We suggest that the discrep- 
ancy is not due to a higher cut-off energy in the 
stellar case, but rather that our analysis is phys- 
ically limited due to a combination of the lack of 
instrument sensitivity, possible power law changes, 
and systematic bias due to temperature effects in 
very small flares (see above). By modeling flares 
as driven dissipative avalanche systems, Lu et al. 
(1993) predict a high-energy rollover in the energy 
distribution whose magnitude depends on the size 
of the active region where the flares occur. Kucera 
et al. (1997) have applied this to SMM/HXRBS 
data and find empirically that the total energy of 
an event has a maximum, Ecuto/j ~ 5 x 10^^^^(jg 
ergs, where ^^^ihs is the total area of active re- 
gions in units of solar micro-hemispheres, i.e., 
l/uhs = 3.04 X 10^^ cm^. Thus, for active stars 
such as the onc^s wc; are considering, where active 
regions covering large fractions of the stellar sur- 
face are expected (A^hg >> 10**), total flare ener- 
gies in excess of 10'^^ ergs are achievable, and thus 
our analyses are valid up to these energies. 

Background Corrections: All the stars in our 

sample are strong sources of EUV emission, and 
the background is generally small compared to 
the source strength. We do not subtract the 
background from the datasets, nor model it sep- 
arately, but assume that the constant component 
in the flare model (Equation 3a) includes contami- 
nation due to the background. That is, we assume 
that the background does not vary on timescales 
smaller than the adopted decay timescale, and 
that any variations that exist in the instrumental 
and astrophysical background are small in magni- 
tude and do not contribute to the flare component. 
Any departures from strict constancy in the back- 
ground will result in a poorer determination of rc 
because of the larger spread in the distribution of 
arrival times, g{dt) (see Figure 5). In cases where 
the above assumptions are invalid, the background 
data will contaminate the signal, but as estimated 
in Tables 1-3, this contamination will be < 10%. 
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3.3.3. Advantages and Disadvantages 

Unlike previous methods that determine flare 
number distributions non-parametrically, that is, 
by directly counting the number of detected flares, 
we model the flare distribution as a power law, 
and generate photon arrival-times to compare with 
the observed event lists. Thus the validity of our 
results is directly dependent on the applicability 
of the model (see extensive discussion above in 
§3.3.2), i.e., the derived parameter values are only 
as good as the model. Thus, our results are correct 
subject to the caveat that we require that all flares 
are physically "similar" and that the flare distri- 
bution follows a power law. But the fact that we 
model the distribution means that we can usefully 
extend the analysis to the regime where the light 
curves are dominated by large numbers of small, 
undetectable flares, and hence this method is best 
suited to study continuous micro-flaring. 

Because the modcil is stochastic in nature, any 
realization of the photon arrival times can be sig- 
nificantly different from the observed event list 
even if the parameter set matches exactly. Thus 
a large number of simulations must be carried out 
for each set of parameters, which is a time con- 
suming process. On the other hand, the stochas- 
tic nature of the model renders unnecessary an ex- 
act match between the model realization and the 
observed light curve, allowing us to explore weak 
flare events. 

As shown in §3.2 (Figure 5) and §3.3.1, the 
method loses sensitivity for a > 3 and loses stabil- 
ity for a < 1.5. However, over the range of a that 
is of scientific interest, i.e., spanning the critical 
value a = 2, the method is robust and can discrim- 
inate between the two cases where flare emission 
si may be a significant contributor to the coro- 
nal emission budget [a > 2), or wlu-re the light 
curve may be dominated by large flares which do 
not contribute significantly to the energy budget 
(a < 2).i3 

'Note that the mere fact that a > 2 does not guarantee 
that the coronal emission is dominated by flaring, but the 
normalization of the power law distribution must also be 
large. Our approach allows us to independently determine 
the flare contribution rp to the total count rate. The re- 
sults (sec §4) show that the normalization is such that the 
flare contribution is generally > 50%. 



4. Results 

We have applied the algorithm detailed above 
(§3) to EUVE/m data on FKAqr, V1054Oph, 
and AD Leo (see §2). The results are summarized 
in Table 4. Below we comment on each dataset in 
detail. 

4.1. FKAqr 

As anticipated above (see Figure 5), the max- 
imum a posteriori value (MAP; see footnote 10) 
for FKAqr are a...„ « 2.5, rp » 0.22 and 
'^'^MAP ^ 0-36, corresponding to a flare contribu- 
tion of K, 65% to the total emission. The joint 
probability distribution of a and r p , marginalized 
over rr, is shown in Figure 7, and the individual 
probability distributions of a and rp, marginal- 
ized over the other parameters, are shown in Fig- 
ure 8. The best-fit values are a = 2.60 ± 0.34, 
rp = 0.19 ± 0.12, and rr = 0.37 ± 0.01. The cor- 
relations between the parameters are evident in 
p{a,rp). Note that there is a small probability 
(-- 0.1) that {rp « rr ^ a < 2}. It is worth 
pointing out the cause of the large values of the 
probability for small rp. First note the presence 
of a secondary peak at [a = 2.2, rp = 0.05) which 
contributes to a larger spread in the uncertainty 
in r^?, and incidentally also indicates the poten- 
tial existence of multiple power-law components. 
Second, this effect is a measure of the stability of 
the model best-fit parameters to the data; if there 
are large parts of the parameter space which pro- 
vide adequate (though not good) fits to the data, 
their contributions are enhanced due to the larger 
volume of the space they occupy. This is indeed 
the case here for combinations of small values oirp 
and a, where the skewness induced in g mod el (St) 
by a < 2.5 is minimized due to the lower weight 
accorded to it because of the smaller values of k. 
Data from Chandra or XMM-Ncwton, with their 
higher expected count rates, are necessary to ex- 
plore the region of smaller values of dt and thus 
better constrain the parameter ranges. 

4.2. V1054Oph 

The light curve of V1054 Oph (Figure 2) shows 
considerable variability with some relatively weak 
flare-like events. This is a signature of a flare dis- 
tribution with large values of a, and indeed de- 
tailed analysis confirms this impression; we find 
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Table 4: Summary of Results 



Star 




Flares Vo" 




P d 
Grange 












logio[ergs s~i] 


[xl028 erg] 


FKAqr 


2.60 ± 0.34 


~ 50 (65) 


2.09 X 10* 


29.48 - 33.54 


9.87 


V1054Opli 


2.74 ±0.35 


70 (85) 


8.6 X 10^ 


29.12 - 32.67 


8.64 


AD Leo 












[2,3] 


2.17 ±0.03 


80 (80) 


2.12 X 10* 


28.99 - 33.66 


8.20 


[4,5] 


2.32 ±0.11 


65 (75) 


2.8 X 10* 


28.99 - 33.52 


3.40 



"Power law index. 

'For the average rate, rW as a percentage of the average total rate; within bra<;kets, for the mode ^F^j^p ■ 
'^Normalization factor for power law distribution, for best- fit parameters (see Equation 6a). 
''Range of flare energies included in calculations (sec Equation 5). 

'^Minimum flare energy to which power law should be extended in order for the flare component to account for the entire emission. 



0.40 F 



FK Aqr: p(a,rp) 




3.5 



0.1 7 r 



a = 2.60±0.34 



0.19 



0.06 - 
0.00 
0.00 



FK Aqr 




2.0 2.5 3.0 3.5 



FK Aqr 







rr=0.18±0.1 1 









0.10 0.20 0.30 0.40 0.50 
r, [ct s"'] 



Fig. 7. — Joint posterior probability distribution 
of a and rp, marginalized over tt, for FKAqr. 
The peak of the distribution lies at a^j^p = 2.5 
and rp = 0.22. 

MAP 

a = 2.74 ± 0.35, with an upper bound that is 
not well constrained (Figure 10). From the joint 
probability distribution p{a,rF), the most prob- 
able values arc a.,.„ ~ 2.5 and rp k, 0.065, 

MAP ^MAP ' 

suggesting that it is likely that almost all the ob- 
served emission originates in flare-like events. 

4.3. AD Leo 

In general, the AD Leo data have been ana- 
lyzed in two separate batches because of the large 
time intervals and large number of counts involved 
(which leads to very long computation times) and 
also because the character of the light curve (see 



Fig. 8. — Marginalized 1-D posterior probability 
distributions of a and rp for FKAqr. Note the 
apparently large contribution at small rp, which 
is a consequence of the fact that the flare contri- 
bution is poorly determined but cannot be ruled 
out for any a for small values of rir. 

Figure 3) appears to change from being domi- 
nated by an intense flare at the beginning (seg- 
ments II+III) to being steadier, with weaker flare 
events (segments IV+V). The light-curves also 
suggest that the flare distribution is character- 
ized by smaller values of a than arc FK Aqr and 
V1054Oph. Detailed analysis confirms the latter, 
with aii+m = 2.17±0.03 and aiv-FV = 2.31±0.11 
(see Figures 11,12) but the distributions p(a) over- 
lap for the two segments, and we cannot rule out at 
the 10% confidence level that the a's are identical 
for the two datasets. However, the trend for the 
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AD Leo (ll + lll) 



a = 2.17±0.03 













1.90 2.00 2.10 2.20 2.30 2.40 



1.00 " 
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0.01 - 

0.00 

0.32 0.34 



AD Leo (ll + lll) 







rF=0.36±0.02 







0.35 0.38 
[ct s'] 



0.40 0.42 



Fig. 9.— As Figure 7, for V1054Oph. 



0.17 ^ 
0.12 - 

0.06 - 
0.00 _ 
1.5 



VI 054 Oph 
a = 2.74±0.35 



2.0 2.5 3.0 



3.5 



0.20 
0.13 

0.07 - 
0.00 . 



VI 054 Oph 



rp=0.05±0.02 



0.00 0.05 0.10 0.15 

r, [ct s"'] 



0.20 



Fig. 10.— As Figure 8, for V1054 Oph. Note that 
there is significant mass in the distribution p{a) 
for large vahics of a, which impHes that the upper- 
bound is not well constrained. However, the lower- 
bound is well determined to be > 2.2. 



analyses of individual segments, (an = 2.03±0.05, 
am = 2.22 ± 0.07, cW = 2.21 ± 0.06, and = 
2.31 ± 0.03) does suggest a gradual steepening of 
the flare distribution when large flares are absent 
from the dataset. 



Fig. 11. — As Figure 8, for segments 2 and 3 of 
AD Leo. The parameters are well determined due 
to the large size of the dataset. Values of a < 2 
can be emphatically ruled out. 



AD Leo 




AD Leo (IV+V) 




0.20 0.25 
[ct s'] 



0.35 



Fig. 12. — As Figure 11, for segments 4 and 5 of 
AD Leo. The best-fit value of a is larger than for 
segments 2 and 3, though the probability distribu- 
tions do overlap. The relatively large spread inrp 
is due to the larger allowed spread in a compared 
to the earlier data. The most probable values are 



= 2.3 and rp. 



0.24, which implies that 



"a 

the data are almost entirely due to flaring emis- 
sion. 



5. Summary 

We have modeled the event arrival times from 

active stars FK Aqr, V1054 Oph, and AD Leo with 
particular attention to the component that arises 
from flare-like events. On the Sun, flares are 
known to be distributed as a power law in energy 
(Hudson 1991; also see Aschwanden et al. 2000 and 



references therein) , and numerous studies have es- 
tablished that strong stellar flares also follow a 
power law distribution, with indices ranging from 
1.5 - 2.5 (see Audard et al. 2000 and references 
therein). This is of considerable interest because 
if the power law index q is > 2, then the coro- 
nal X-ray losses could in principle be ascribed to 
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weak flare events that are nevertheless numerous 
enough to dominate the emission. 

We consider active dMe stars with known vari- 
abihty in their light curves where numerous flares 
are seen. Audard et al. (2000) find that in gen- 
eral dF and dG stars tend to have a > 2 while 
dK and dM stars tend to have lower a. In par- 
ticular, they analyze an older dataset of AD Leo 
(from May 1996) and show that detectable flares 
have a e [1.18, 2.35] or a = 2.02±0.28 using differ- 
ent methods. A more detailed analysis by Giidel et 
al. (2002) based on a larger sample of the dataset 
used here shows 2.0 < a < 2.5 . 

We model the event arrival times using a sim- 
ple two-component model comprising of a con- 
stant rate component and a statistical ensemble 
of flare components, with the flare energies dis- 
tributed as a power law. In general, the simplify- 
ing assumptions we make (e.g., constancy of de- 
cay timescales, ignoring the rise times, assuming 
a constant counts-to-energy conversion factor, in- 
cluding the background directly in the model, etc.) 
arc conservative, and tend to underestimate the 
value of a. We find that all the stars in our sam- 
ple clearly have a > 2: for AD Leo, a lies in the 
range 2.06 - 2.32; for FKAqr, a = 2.60 ± 0.34; 
and for V1054Oph, a = 2.74 ±0.35. We thus con- 
clude that coronal heating on these stars is dom- 
inated by impulsive energy release events whose 
energy output is > 2 — 3 x 10^^ ergs, reaching to 
the microflare range. These results are in contrast 
to the solar case, where over similar flare energy 
ranges the observed distribution of flares is shal- 
lower, with a w 1.8, i.e., below the critical value 
of 2. 

Further, we directly estimate the contribution 
of the flare emission to the total observed count 
rate, as one of the parameters deflning the model. 
Because the energy range over which the model is 
defined spans over 4 orders of magnitude, and the 
power law indices indicate steep distributions, we 
expect that the flare component should contribute 
significantly to the total emission. Indeed, we find 
this to be generally > 50%, and in some cases 
being > 80%. Thus, there appears to be no truly 
"quiescent" emission on some of these low-mass 
active stars, i.e., emission from apparently stable 
active region loops as on the Sun is not a dominant 
component of the observed emission. 

We have also explored the possible dependence 



of the various model parameters on flare energies. 
The long observation of AD Leo suggests that a 
increases when strong flares are not evident in the 
data suggesting that the flare distribution steepens 
for flares of smaller energies, though we cannot 
rule out the possibility of a statistical fluctuation 
that mimics this trend. We have also searched for, 
but do not find, evidence of strong dependence of 
the decay timescale on flare energies. 

We flnd that if the flare distributions extend 
to the microflare regime, the energy output due 
to these weak flares is quite sufficient to account 
for the entire coronal emission in the EUVE/DS 
passband. However, it must be noted that the er- 
ror bars on the parameters derived for the fainter 
stars FKAqr and V1054Oph are quite large (for 
instance, values of a < 2 cannot be completely 
ruled out for FKAqr, and a firm upper bound 
on a cannot be set for V1054Oph) and it would 
be of considerable interest to verify and improve 
these results (and also to extend them to a larger 
sample of stellar types) using high-quality data 
such as those obtainable with Chandra and XMM- 
Newton. Data from these observatories are char- 
acterized by good time resolution and in general 
larger count rates, and will therefore allow us 
to explore the arrival time difference distribution 
functions g{5t) at smaller values of 6t, thereby ex- 
tending the range of values of a that the method 
is sensitive to. 
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A. Glossary of terms 

Here we compile, for reference, all the symbols used in the text (Table 5). 



Table 5: Glossary of terms 



Symbol Description 



First Use 



r(x) 

AT 

e(.) 

a 

5t 

K 

m 
pi 

T 

x' 

C{t) 
D 
E 
Fj 

M,M' 

Nf 
I 

dN 
c 

^max 

fit) 

gm 

k 

p{A\B) 

r{t) 

rc 

tf 

rr 

var{x) 

X 

X \ rt 



the Gamma function, whicli for integer a; is (a; — 1)! 
total duration of observation 
Heaviside step function 

power- law index 

arrival time difference between two consecutive photons 
normalization of the power-law in counts units 
Instrument correction factor that includes Primbsching, etc. 
the fraction of time a source spends at rate r, 
flare decay time scale 

= 5^ ( °^rror{pata\' ) ' ^ Statistic measuring the quality of a fit 

model light curve, observable counts in [f, t -{■ dX\ 

observed data 

energy output of flare event 

peak intensity of model flare j 

set of parameters defining the model 

total number of fiares in the model 

information necessary to define the problem 

number of flare events in [£, E + dE\ 

counts due to a flare 

maximum model counts possible due to a flare, for a given dataset 
minimum model counts due to a flare that is practicable to consider 
flare model light curve intensity at time t 

frequency histogram of the distribution of &t 
normalization of the power-law in energy units 

probability that statement A is true given that statement B is true 

model light curve intensity 

constant component model count rate 

mean model count rate of flare component over duration of observation 
mean total model count rate over duration of observation 
variance of quantity x 
mean value of quantity x 

maximum a posteriori value of x, where p{x) is maximum 



§3.2, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§1, Eqn. 1 
§3.2 

§3.1, Eqn. 
§3.1, Eqn. 
§3.2, Eqn. 
§3.1, Eqn. 



10 
3c 
3b 



4 

3a 
9b 

3b 



§3.2, Figure 5 



§3.1, Eqn. 
§3.2, Eqn. 
§1, Eqn. 1 
§3.1, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.2, Eqn. 
§1, Eqn. 1 
§3.1, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.2, Eqn. 
§3.1, Eqn. 
§3.2, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.1, Eqn. 
§3.2, Eqn. 
§3.2, Eqn. 



3a 
7a 

3b 

2 

3b 
7a 

4 
5 
5 

3b 

9a 

6a 

7a 

3a 

2a 

2a 

2b 

12a 

12a 



§3.2, Footnote 10 
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